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Abstract 

In this work, we revisit the question of stability of multibreather configurations, i.e., discrete 
breathers with multiple excited sites at the anti-continuum limit of uncoupled oscillators. We 
present two methods that yield quantitative predictions about the Floquet multipliers of the 
linear stability analysis around such exponentially localized in space, time-periodic orbits, based 
on the Aubry band method and the MacKay effective Hamiltonian method and prove that their 
conclusions are equivalent. Subsequently, we showcase the usefulness of the methods by a 
series of case examples including one-dimensional multi-breathers, and two-dimensional vortex 
breathers in the case of a lattice of linearly coupled oscillators with the Morse potential and in 
that of the discrete 4 model. 

1 Introduction 

Over the past two decades, there has been an explosion of interest towards the study of Intrinsic 
Localized Modes (ILMs), otherwise termed discrete breathers [1]. This activity has been, to a 
considerable extent, fueled by the ever-expanding applicability of these exponentially localized in 
space and periodic in time modes. A partial list of the relevant applications includes their emergence 
in halide-bridged transition metal complexes as e.g. in [2], their potential role in the formation of 
denaturation bubbles in the DNA double strand dynamics summarized e.g. in [3], their observation 
in driven micromechanical cantilever arrays as shown in [4], their investigation in coupled torsion 
pendula [5], electrical transmission lines [6], layered antiferromagnetic samples such as those of a 
(C2H5NH3)2CuCi4 [7], as well in nonlinear optics [8] and possibly in atomic physics of Bose-Einstein 
condensates [9, 10] and most recently even in granular crystals [11]. 

In parallel to the above experimental developments in this diverse set of areas, there has been a 
considerable progress towards the theoretical understanding of the existence and stability properties 



of such localized modes summarized in a number of reviews and books; see e.g. [1, 8, 12, 13]. 
Arguably, one of the most important developments in establishing the fundamental relevance of 
this area in coupled nonlinear oscillator chains has been the work of MacKay and Aubry [14], 
which established the fact that if a single oscillator has a periodic orbit (and relevant non-resonance 
conditions are satisfied), then upon inclusion of a non- vanishing coupling between adjacent such 
oscillators, an ILM type waveform will generically persist. 

Given the confirmation of persistence of such modes, naturally, the next question concerns their 
robustness under the dynamical evolution of the relevant systems, which is critical towards their 
experimental observability. This proved to be a substantially more difficult question to answer 
in a quantitative fashion, especially so for ILMs featuring multiple localized peaks, i.e., multi-site 
breathers (since single-site breathers are typically stable in chains of linearly coupled anharmonic 
oscillators). Two principal theories were proposed for addressing the stability of such periodic orbit, 
discrete breather states (and identifying their corresponding Floquet multipliers). Interestingly, 
these originated independently from the same pioneers which established (jointly) the existence of 
such modes in [14]. In particular, the first theory was pioneered by Aubry in his seminal work of 
[12] and will go under the name Aubry Band (AB) theory, hereafter. The second one is an effective 
Hamiltonian method which was introduced in a series of papers by MacKay and collaborators [15] 
(and will be termed accordingly MacKay Effective Hamiltonian method (MEH)). The AB approach 
was adapted to the stability of discrete breathers and multibreathers in the setting of Klein-Gordon 
lattices in the work of [17]; see also [18]. The MEH approach was applied to the same setting in 
the recent work of [19]; see also [20]. 

Our aim in the present work is to unify the two methods by firmly establishing the equivalence 
of the stability conclusions of the Aubry band and MacKay effective Hamiltonian methods. Sub- 
sequently, we illustrate the usefulness and versatility of the methods, we apply them to a range of 
physically interesting chains of oscillator model examples, such as the Morse potential which arises 
in the study of DNA bubbles [3], as well as the </> 4 potential which arises in applications in dusty 
plasmas [21], as well as in field theory, particle physics and elsewhere; see e.g. the recent discussion 
of [22] and the earlier review [23] and references therein. Our presentation is structured as follows. 
In section 2, we compare the two approaches and showcase the equivalence of their conclusions. 
In section 3, we study multibreathers and vortices in the case of the Morse potential. In section 
4, we present the corresponding results for the 4> 4 Klein-Gordon lattice. Finally, in section 5, we 
summarize our findings and present our conclusions. 



2 Comparison between the two approaches 
2.1 Preliminaries - Terminology 

The relevant system under consideration will be a Klein-Gordon chain of oscillators with nearest- 
neighbor interaction and Hamiltonian 



H = Hr, + eHi 



oo 
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',= — oo 



pi + v{ Xi 



oo 
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(1) 



As indicated previously, we will examine the two approaches (AB and MEH) for the linear 
stability of multi-site breathers of this general class of systems. Both approaches are based on the 
notion of the anti-continuum limit. In this limit (e = 0) we consider n "central" oscillators moving 
in periodic orbits with the same frequency oo (this will be our "multibreather" for e / ), while 
the rest lie at the equilibrium (x,x) = (0,0). For e ^ some of these configurations, depending 
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on the phase differences between the oscillators, are continued in order to provide multibreather 
solutions. It is interesting/relevant to note here that while the MEH approach provides explicit 
conditions about which configurations can be continued to finite e (the critical points of the relevant 
effective Hamiltonian) , the AB theory provides only stability information for a given configuration 
(for which we already know otherwise that it should exist at finite e). 

The linear stability of these solutions is determined by the corresponding Floquet multipliers. 
For a stable multibreather we require that all the multipliers lie on the unit circle. In the anti- 
continuum limit these multipliers lie in three bundles. The two conjugate ones, that correspond to 
the non-central oscillators, lie at e izlU)T , while the third one lie at +1 and consists of n multiplier 
pairs, corresponding to the central oscillators. Each pair of +l's correspond to the phase mode and 
growth mode of each isolated excited oscillator, meaning that a small change in the initial phase 
or a small change in frequency leads to an extremely close periodic solution, for the growth mode 
with slightly larger or smaller amplitude. 

For e 7^ 0, the non-central corresponding bundles split and their multipliers move along the unit 
circle to form the phonon band, while the multipliers at unity can move along, either the unit circle 
(stability), or along the real axis (instability). However, a pair of multipliers always continues at 
+1 corresponding to the phase mode and growth mode of the whole system. Hence, the stability of 
the multibreather, at least for small values of the coupling, is determined by the multipliers of the 
central oscillators. For larger values of e, a Hamiltonian-Hopf bifurcation can occur and destabilize 
an initially stable multibreather. 

At this point, it is relevant to make a note in passing about the striking similarities between 
the discussion above (at and near the anti-continuum limit) with that of the linear stability of 
standing waves in the discrete nonlinear Schrodinger (DNLS) equation. In that case, due to the 
monochromatic nature of the solutions and the U(l) invariance of the latter model, it is possible 
to directly consider the eigenvalues associated with the standing wave solutions. However, there 
is a direct analogy with the spectrum of the excited sites being associated with the eigenvalues 
at the origin at the anti-continuum limit and the continuous spectrum lying at a finite distance 
from the spectral plane origin, and how at finite coupling these zero eigenvalue pairs of the excited 
oscillators are the ones that may give rise to instability. In fact, it turns out that even the conditions 
under which instability will ensue for multibreathers of the KG directly parallel the ones for multi- 
breathers (or multi-site standing waves) of the DNLS. The latter are analyzed in considerable detail 
for 1—, 2— and 3— dimensional settings in [13]. 

Returning to our KG setting, the MEH approach considers the Floquet multipliers given as 
A = exp(crT), with T = 2n/uj, whereas in the AB approach, A = exp(i6>). Then, 

Due to the symplectic character of the Floquet matrix if A is a multiplier, so is A -1 , and due 
to its real character if a A is non real multiplier, so is A*, where the asterisk denotes the complex 
conjugate. Therefore, the corresponding multipliers come in complex quadruplets (A, A -1 , A*, A* -1 ) 
if |A| / 1 and A is not real, or in duplets A, A -1 if A is real, or A, A* if |A| = 1 and not real. In 
addition, due to the time translation invariance of the system there is always a pair of eigenvalues 
at +1. This has as a result that both cij's and #j's, come also in quadruplets or duplets (a, —a) 
if a is real, ((6*,— 6*) if 9 is imaginary) or (a*, — a*) if a is imaginary ((9,-6) if 9 is real). In 
principle, the duplets could collapse at a single value A = ±1, but there is always a pair of +l's for 
the systems under study, as explained above. 
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2.2 The MEH approach 



The MEH approach consists of constructing an effective Hamiltonian, whose critical points are 
in correspondence with periodic orbits (in our case multibreathers) of the original system. This 
method has been originally proposed in [15] and used in the present form in [16]. The effective 
Hamiltonian can be constructed as follows. 

After considering the central oscillators we apply the action-angle canonical transformation to 
them. Note that, in the anticontinuous limit, the motion of the central oscillators, in the action- 
angle variables, is described by Wi = corf + Wio, Ji = const., for i = 0, ... , n — 1. Where Wi is 
the angle, Wio is the initial phase and Ji the action of the i-th. central oscillator. For this kind of 
systems, the action of an oscillator can be calculated as 

Jl = 7 LJ PidXi = 7 kJ foWl 2 * 1 *- ( 3 ) 

Since we are interested in a first order approach, the effective Hamiltonian can be written as 
H cS = Hq + e(Hi), by neglecting terms which do not contribute to the results in this order of 
approximation. In this formula, (H\) is the average value of the coupling term of the Hamiltonian, 
over an angle in the anticontinuous limit, which is equivalent to the average value of H\ over a 
period 

(Hi) = ^ I H^t. 



This averaging procedure is performed in order to lift the phase degeneracy of the system. For the 
same reason we introduce a second canonical transformation 

•& = w A = Jo + • • • + J n -i 

n—l 

4>i = Wi- Wi-i Ii = ^2jj i = 1, . . . , n - 1. 

j=i 

In these variables, the effective Hamiltonian reads 

H cS = H (I i ) + e(H 1 )(4> i ,I i ). (5) 

Note that, since the calculations are performed in the anticontinuous limit, the contribution of the 
non-central oscillators has disappeared. 

As we seldom know the explicit form of the transformation (x,p) i->- (w,J), we use the fact 
that since the motion of the central oscillators for e = is periodic, and possesses the t h-> — t, x i— > 
x,p I—?- —p symmetry, it can be described by a cosine Fourier series Xi{t) = XlfcLo Ak(Ji) cos(kwi). 

Note that at the anti-continuous limit, the orbits differ only in phase (i.e. Wj = oj Vi), therefore 
Ji = J and the coefficients A^s do not depend on the index i. 

So, excluding the constant terms, (Hi) becomes for the KG problem [19] 

oo n—l 

(Hi) = --J2J2 A l co < k ^ ( 6 ) 

k=l i=l 

One of the main features of the MEH approach is that the critical points of this effective 
Hamiltonian correspond to the multibreather solutions of the system. This fact provides the cor- 
responding persistence conditions, as the simple roots of 9( £i^ = 0. Remarkably, in this setting, 
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similarly to what is known also for the DNLS [13], it can be proved that the only available multi- 
breather solutions in the one-dimensional case are the ones with relative phase among the excited 
sites of or it. 

The second important fact the MEH approach yields is that the linear stability of these critical 
points (i.e., the Hessian of the effective Hamiltonian) determines the stability of the corresponding 
multibreather. In particular, the nonzero characteristic exponents of the central oscillators u\ 
(see the discussion in the previous subsection) are given as eigenvalues of the stability matrix 

-I 

1 O 



E = JD 2 H eS where J = 
in (5) for the H eS we get: 



is the matrix of the symplectic structure. By using the form 
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Since the only permitted values of the relative phases are fa = 0, or fa = tt, the matrix simplifies 
considerably acquiring the form: 



E = 



O B 
C O 



O eBi 
C + eCi O 



(8) 



which, subsequently, if we consider only the dominant eigenvalue contributions, we get that of 
eX-BC; where xbc ar e the eigenvalues of the (n-lxn-1) matrix Bi • Co which reads 



du dco 

Bi • C^n = ^ — Z = — 

dJ dJ 



( Vi -h 
-h 2/ 2 



V 





-h o 

-fn-2 2/ n _ 2 
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\ 



— fn-2 
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(9) 



In this expression oj = dHo/dJ denotes the frequency, while 

1 oo 

fi = m) = -^2 k2A l cos ( k ^)- 



(10) 



k=l 



This leads to the characteristic exponents (i.e., effective eigenvalues) of the DB in the form: 



with Xz being the eigenvalues of the (n — 1 x n — 1) matrix 



(11) 



r Zi 



-fi 



Z%^% — 2/j 
otherwise. 



(12) 
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2.3 The AB approach 

We demonstrate hereby that (11) can be reobtained based on the AB approach, by using the 
exposition of [17]. To this end, we recall that the aim of the AB approach is to look for the 
displacement that Aubry's bands [12] experience when the coupling e is switched on. What we 
plan to do below is to calculate the Floquet eigenvalues assuming that the bands are parabolic and 
their shape does not change when the coupling is introduced. 

First, we recall the basics of Aubry's band theory with the notation used in [17] adapted to the 
notation in the present paper, where convenient, for ease of comparison. The Hamilton equations 
applied to the Hamiltonian of Eq. (1) can be written as: 

x n + V'[xi) + e-^- = i = l,...,N, (13) 
for a generic coupling potential Hi, or, if it is harmonic: 

N 

x n + V'(xi) + e C ijXj = i = l,...,N (14) 
i=l 

where C is a coupling constant matrix. Let us define x = [xi(t), . . . , XN(t)]^ (t meaning the 
transpose matrix). Defining V(x) = [V(xi), . . . , V(xn)]\ dHi/dx = [dHi/dxi, . . . , OHi/Oxn]^ 
and so on, Eq. (13) can be written as: 

x + V\x) + e^- = 0. (15) 

Suppose that x(t) is a time-periodic solution, with period T and frequency oj, its (linear) stability 
depends on the characteristic equation for the Newton operator M e given by 

A4(uK ='(+ V"{x)H + e^£= Et, (16) 

where * product is the list product, i.e., f(x) * £ is the column matrix with elements f(xi(t)) ^«(t), 
and d 2 Hi/dx 2 is the matrix of functions d 2 Hi/dxidxj, which depends on t through x = x(t). 

If E = 0, this equation describes the evolution of small perturbations £ = £(t) of x = x(t), 
which determines the stability or instability of x. It is however, extremely useful to consider the 
characteristic equation for any eigenvalue E as it is the cornerstone for Aubry's band theory. 

Any solution £ of Eq. (16) is determined by the column matrix of the initial conditions for 
positions and momenta Q(0) = [£i(0), . . . £jv(0),7ri(0), . . . 7Tjv(0)]t, with 7Tj(t) = £i(t). A basis of 
solutions is given by the 2 iV functions with initial conditions 1F(0), v = 1,...,2N, with 1^(0) = 

The Newton operator depends on the T-periodic solution x(t), and therefore, it is also T- 
periodic and its eigenfunctions can be chosen also as eigenfunctions of the operator of translation 
in time (a period T). They are the Bloch functions i{6i,t) = x(&i,t) exp(i#j t/T), x{@i~k) being a 
column matrix of T-periodic functions. The sets {£(#«, 0), £(#«, 0)} are also the eigenvectors of the 
Floquet operator Te or monodromy, that maps f2(0) into 0,(T), that is, fJ(T) = Fe£1(0). Their 
corresponding eigenvalues are the 2N multipliers {Aj} = exp(0j), with {9i} being the 2N Floquet 
arguments. 

The set of points (6,E), with 6 being a real Floquet argument of Te, has a band structure. 
As the Newton and Floquet operators are real, the Floquet multipliers come in complex conjugate 
pairs. Therefore, if (9,E) belongs to a band (i.e. is real), {—9,E) does it too, i.e., the bands are 
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symmetric with respect to 8, which implies that dE/d8(0) = 0. There are always two T-periodic 
solutions, with Floquet multiplier A = 1 (8 = 0) for E = 0. One is x(t), which represents a change 
in phase of the solution x(t) and it is called the phase mode; the other is called the growth mode, 
given by dx(t)/dui, and represents a change in frequency and consequently in amplitude. The 
consequence is that there is always a symmetric band tangent to the axis E = at 9 = 0. 

There are at most 2N points for a given value of E and, therefore, there are at most 2 N bands 
crossing the horizontal axes E = in the space of coordinates (8,E). The condition for linear 
stability of x(t) is equivalent to the existence of 2N bands crossing the axis E = (including 
tangent points with their multiplicity). If a parameter like the coupling e changes, the bands evolve 
continuously, and they can lose crossing points with E = 0, leading to an instability of the system. 

The first item to find out are the bands at the anticontinuous limit, where Eq. (14) reduces to 
N identical equations: 

x l + V'{x t ) = 0. (17) 

If we consider solutions around a minimum of V, the oscillators can be at rest Xi = 0, or oscillating 
with period T; the latter are identical except for a change in the initial phase, so they can be written 
as Xi(t) = g(cot + Wio) with g(ut) being the only T-periodic, time-symmetric solution of Eq. (17) 
with g(0) > g(7r). Therefore, the excited oscillators can be written as: 

OO OO OO 

Xi(t) = zq + 2 ^2 z k cos[k(ujt + Wio)] = ^2 Ak cos[k(ut + Wio)] = ^ cos(kwi) , (18) 

k=l k=0 k=0 

with Ak = 2zk if k > 0, Aq = zq and Wi = cot + wiq. 

Let n be the number of excited oscillators at the anticontinuous limit, labeled i = 0, . . . , n — 1. 
Then, there are n identical bands tangent to the axis E = at 8 = for each excited oscillator, 
and N — n bands, corresponding to the oscillators at rest, with 2(iV — n) points intersecting the 
E = axis. 

Thus, the excited bands can be approximated around (8, E) = (0, 0) by 

E{6) ^E + k8 2 , (19) 
with Eq = e\ q and Xg being the eigenvalues of the (n x n) Q-matrix defined below. Additionally, 

1 d 2 E oj 2 dH 1 dH 

where we have made use of [17, Eq. (B14)]. The factor k is positive if the on-site potential V is 
hard and negative if V is soft (a potential is hard if the oscillation amplitude increases with the 
frequency and soft otherwise). When the coupling is switched on, the bands will move and change 
shape; the E = eigenvalue is degenerate with multiplicity iV — n at e = 0, but this degeneracy 
is generically lifted for e / and only one band will continue being tangent at (8, E) = (0, 0) due 
to the phase mode. Applying degenerate perturbation theory to Eq. (16), with eH\ being the 
perturbation, a perturbation matrix Q can be constructed [17], whose eigenvalues Xq are those of 
the perturbed Newton operator. The non-diagonal elements of Q are given by 

1 f T d 2 H 

Qij = / ii- — ^-xjdt, i^j, i = 0...n-l, j' = 0...n-l, (21) 

Hi Hj .A, dxi dxj 



with ^ = \J \xi) 2 dt. The diagonal elements are 



Qn = -Y J -Qi3- (22) 



7 



If the on-site potential V(xi) is homogeneous and the coupling is given as in Eq. (1), as is the case 
in the present paper, jUj = (2ttJ) 1 ^ 2 Mi. Let us calculate the derivatives of Hi, hij = d 2 H\/dxi dxj. 
Because of the way the diagonal elements of Q are constructed, we only need the derivatives with 
i / j. It is easy to see that they are zero except for hi-ij = h^-i = —qi (defined below) for 
i = 1, . . . ,n — 1. The derivatives /io,n-i and /i n -i,o are also zero as the oscillators at the extremes 
of the multibreather are not coupled between them. Then, the matrix Q becomes: 

Qi,i-\ = Qi-i,i = -Qi, fori = l...n-l 
Qi,i = 5i-i + for i = 1 ... n - 2 

Qo,o = 5i (23) 

Qn—l,n—l — 5n— 1 

otherwise 
or, explicitly: 



Qi,j — < 



/ qi -qi 
-?i qi + 52 -52 



Q 



~q n -2 5n-2 + 5n-l ~5n-l 

-5n-l 5n-l / 



(24) 



with 



g, ee g (&) = ^ = £ E k 2 A\ cos(^) = » fi , * = l,...n-l,, (25) 



/ J [ii(t)] 2 dt 2 J 



k>i 



Then, by using [26, Lemma 5.4] we see that the matrices Q and — Z have the same nonzero 



eigenvalues i.e. 



Xq = ~j%-Z 



In addition, Q has also a zero eigenvalue. 

Some important values of q((f>) are the following ones: 



5(0) = 1 



E fc >i(-i) fc fc% 2 

Sfc>l ^ 2 ^fe 



= - 7 



(26) 

(27) 
(28) 



For a Morse potential, 7 = w; for an even potential, 7 = 1. 

According to the AB theory [12], the Floquet multipliers are given by the cuts of the bands 
with the E = axis; thus 



K 



kJ 



Xz 



and, applying the last results 



= ±T x! ep jXz , 



(29) 
(30) 
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where we have taken into account that 



i 



dH _dHdJ _ dJ 
duo dJ du du 

Finally, introducing (30) into (2), we get (11), which completes the proof of equivalence of the 
relevant Floquet multiplier predictions. 



3 The Case Example of the Morse potential 
3.1 Characteristic exponents 

We now consider some special case examples, starting with a linearly coupled lattice of oscillators 
subject to the Morse potential. As indicated previously, the only configurations that may exist 
in the one-dimensional setting are ones which involve excited oscillators either in-phase (i.e., with 
4>i = 0) or out-of-phase (i.e., with = 7r); see [18], [19] and also [24] for a detailed discussion. 
Here, we proceed to perform some explicit calculations for the Floquet multipliers a in the case 
of n-site breathers. In what follows we consider only the positive a. To this end we express (11) 
making use of (26) : 2 

* = ftt*v> (32> 

where Xq(4>) denotes the Q-matrix eigenvalues for a given 4>. It is straightforward to show that 3 

X 9 (0) = 4sin 2 ^ m = l,...,n-l (33) 

and that 



X,M = -47X 9 (0) (34) 

For instance, in the case of a 2-site breather, x<j(0) = 2 and x<?( 7r ) = —^l- 

We now focus on the particular case of the Morse potential, since it is a potential for which 

closed form analytical expressions can be found. [For other types of potentials, some approximations 

can be made for small and high frequencies; alternatively, the required single-oscillator parameters, 

such as J and du/dJ can be calculated numerically]. 

In the Morse case and in order to evaluate J and dco/dJ, we express J as a function of the 

Fourier coefficients: 



J = 2uj^2k 2 zl (35) 



k>l 



For this potential, 



lr The expression w = dH/dJ comes from the Hamilton's equations for the action-angle variables, since all the 
calculations are performed in the uncoupled, and therefore integrable, limit. 

2 In what follows, and in order to fix ideas, given the equivalence of the two methods, we will use the formulation 
with the Q-matrix. 

3 We are neglecting the eigenvalue, associated with m = 0. 
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Substituting into the action 



J = l-w->|£ = -l. (37) 



Thus, 



cr(0 = \ e 



1-w 



Co' 



X.WO- (38) 



In the case of a general phase, we can express Xg(0) = g(</>)x<?(0) with c/(0) given by (25). In 
the special case of the Morse potential, we have 

QW = ^5> fe cos(£^). ( 39 ) 

k>l 

To obtain the relevant sum, we use a simple geometric series formula that can be found e.g. in 
[25], according to which: 

/ ±\ 2w cos^-r 

g(0) = — r ^. (40) 

J 1 - 2rcos^ + r 2 V ' 

Consequently, 



^) = V^ i-"«.V + ^ (0) - (41> 



For the relevant values of </> for time-reversible multibreathers, we get: 



, 1 — U> , . m7T / 1 — LJ ... 

^(0) = ye-^X?(0) = 2sm — ^e—^ m = l,...,n-l (42) 

<t(7t) = ^/-e(l - w)x,(0) = 2sin^V-e(l-^) m = l,...,n-l (43) 

Figures 1 and 2 show, respectively, the analytical eigenvalue predictions (dashed lines) for stable 
and unstable two-site and three-site breathers with the Morse potential and how they favorably 
compare to the corresponding numerical results (solid lines), obtained via a fully numerical linear 
stability analysis (and corresponding computation of the Floquet multipliers). It is clear that the 
predictions are very accurate close to the anti-continuum limit, and their validity becomes progres- 
sively limited for larger values of the coupling parameter e, yet they yield a powerful qualitative 
and even quantitative (in the appropriate parametric regime) tool for tracking the stability of these 
localized modes. The figures also illustrate typical profiles of the corresponding two- and three-site 
ILMs. 



3.2 Vortices in square lattices 

The methodology can also be extended to lattices of higher dimensionality. We consider below some 
basic properties of discrete vortex breathers of different integer topological charges, in a square 2D 
lattice. Firstly, we consider square vortices over a single "plaquette" of the 2D lattice with 5 = 1, 
i.e., at the anti-continuum limit, the excited sites are (0,0), (0,1), (1,1) and (1,0) with a phase 
difference <f> = tt/2 between nearest neighbors. This implies a perturbation matrix given by: 
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0.01 0.02 0.03 0.04 0.05 0.01 0.02 0.03 0.04 0.05 



e e 

Figure 1: (Top panels) Profiles of an in phase (left) and an out-of-phase (right) 2-site breather with 
the Morse potential for to = 0.8 and e = 0.05. The bottom panels show the value of the characteristic 
exponents a of the corresponding configurations, with respect to the coupling parameter e. Dashed 
lines correspond to the predictions of the stability theorems, while solid ones to full numerical linear 
stability analysis results. 



Q = g(7r/2)Q,with Q 



(2 -1 -1 \ 
-12-10 
0-12-1 
\ -1 -1 2 J 

with g(-7r/2) given from (40) which is evaluated as: 



q(n/2) 



J(l+u; 2 ) 

This, in turn, upon use of Eq. (32) implies that 



a = i(l — uS) 



(44) 



(45) 



(46) 



1 + w 2 " 

where Xa are the eigenvalues of the Q matrix. This corresponds to the matrix of the normal 
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Figure 2: Same as in Fig. 1, but now for the unstable left-panel configuration of three in-phase 
excited sites (with two real multiplier pairs as shown in the bottom panel) and the potentially stable, 
close to the anti-continuum limit, case of the out-of-phase, three-site right-panel configuration. 



modes of a ID chain of 4 linearly coupled oscillators with periodic boundary conditions. Let us 
recall that for a system of n coupled oscillators, the eigenvalues are given by: 



■ 2 m7r , 

Xo=4sm m=l,...,n — 1, (47) 

n 

in addition to the eigenvalue. In the present case of 4 oscillators with periodic boundary con- 
ditions, the nonzero eigenvalues are given by 2 and 4, with the former being doubly degenerate. 
Thus, we have for 5 = 1 vortices the following spectrum: 



i(l — uj)y2 single eigenvalue 

(48) 

2i(l — u) \J j^ji double eigenvalue 

which implies stability for e > 0. 

This type of analysis can be generalized for arbitrary values of the vorticity S, leading to the 
conclusion that Q is the matrix of 45 coupled oscillators, which implies that vortices with any 
integer topological charge will be stable for e > in the case of a lattice with an on-site Morse 
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potential. For instance, in the case of the S = 2 vortex, we obtain the explicit expressions for the 
eigenvalues: 



a = < 



' ^1 /(2-2 1 / 2 ) 



1+UJ 2 



single eigenvalue 



i(l-u>) 
2 



j^jz double eigenvalue 



(49) 



i(l-w)^ 



double eigenvalue 



single eigenvalue 



It is important to highlight here some interesting differences between the above results and the 
case of the DNLS (and more generally that of even potentials in KG chains, including the case of 
the hard 4> 4 lattice considered below). In the latter class of problems, the vanishing of the odd 
coefficients in the Fourier expansion of the periodic orbit leads to the conclusion that q(ir/2) = and 
hence there is no contribution to the eigenvalues to leading order. This is the situation which has 
been characterized as "super-symmetric" in [13] and one in which the higher order contributions 
would be critical in determining the stability. Nevertheless, in the case considered herein, the 
asymmetry of the Morse potential produces a nonvanishing of q(n/2) and offers a corresponding 
nonzero leading order correction to the eigenvalues at 0(e 1//2 ). 

Figure 3 shows the dependence of stability eigenvalues for the S = 1 and S = 2 vortices and their 
comparison with the obtained fully numerical linear stability results as a function of the coupling 
e. As can be observed in the figures, the approximation is less accurate in this case, although it 
is qualitatively correct. The reason for the partial disparity is that higher order contributions to 
the relevant eigenvalues (whose calculation is considerably more technically involved) lead to the 
observed splitting of all the doubly degenerate eigenvalue pairs. In the relevant cases, the analytical 
(dashed line) predictions can be seen to straddle the two observed numerical pairs. 



4 The Case Example of the Hard 4 Potential 

The time evolution of a single oscillator in the hard (f) 4 potential, V(x) = x 2 /2 + x 4 /4 is given by: 



2m f t \ 2m (2K(m) 
X W = \ cn I =, m = \ cn I ut, m , (50) 



I -2m VV1 - 2m' / V 1 - 2m \ vr 
where cn is a Jacobi elliptic function of modulus m and K{m) is the complete elliptic integral of 
the first kind defined as K(m) = f^ 2 [1 — msin 2 x] -1 / 2 dx. 

The breather frequency u is related to the modulus m through: 

LJ= ; = ^— -■ (51) 

2y/l - 2mK(m) v ; 

The elliptic function can be expanded into a Fourier series leading to [28]: 



TT 



Z ^ 1 = K(^)\lT^r^^ , = 0,1,2,.... (52) 



where q is the elliptic Nome which is defined as 
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Figure 3: The characteristic exponents of vortex configurations with 5 = 1 (left) and S = 2 (right), 
with respect to e, for the Morse potential and uj = 0.8. Dashed lines correspond to the theoretical 
predictions based on Eqs. (48) and (49), respectively; the full numerical linear stability results are 
given by solid lines and indicate that all doubly degenerate eigenvalue pairs split due to higher 
order contributions in the relevant expansions in the coupling constant e. 



q = q(m) = exp(— -kK(1 — m)/K{m)). (53) 

In order to get X 9 (0) and <9w/<9J, we cannot use (10) and (25) because it is not possible to find 
a closed form expression. Instead, we use the integral expression: 



1 



r 



/(&) = IT- / Xi(t)x i+l (t)dt. (54) 
After some manipulations (where it is crucial to apply [27, identity 171]), we obtain: 



8K(m) 



7r 3 a;(l — 2m) 

8K(m) 
7r 3 a;(l — 2m) 



[cs(a, m)ns(a, m)[2E(m) — K(m)(l + dn 2 (a, m) 
[ds(a, m)(cs 2 (a, m) + ns 2 (a, m))Z(a, m)] 



(55) 



where E(m) is the complete elliptic integral of the second kind defined as E(m) = j^ 2 [1 



msin 2 x] 1 / 2 dx, Z(a,m) is the Jacobi zeta function and a = 2K(m 
For the action J, a similar manipulation leads to 



7T. 



J 



16K(m) 
3vr 2 



1 



m 



1 - 2m 



K(m) — E(m) 



(56) 



The derivative of this expression is cumbersome to handle. So, in what follows, we will work 
instead with numerically obtained values of J and dui/dJ which are relevant for time-reversible 
multibreathers and vortex breathers, as for these cases we need /(0), /(tt) and f(n/2). As indicated 
previously, for every even potential, Z2 U +i = 0, and, consequently, /(0) = J/u, f(n) = — /(0) and 
/(vr/2) = 0. This leads to: 
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Figure 4: Dependence with respect to u) of the action (left) and dto/dJ (right) for the hard <\> 
potential. The dashed line corresponds to the prediction of the RWA [Eq. (59)], while the solid 
one represents the exact numerical result. 



Figure 4 shows the dependence of J and dto/dJ with respect to the frequency. Figures 5 and 6 
illustrate subsequently the relevant stability eigenvalues for 2-site and 3-site breathers as obtained 
from the expressions above and compare them to the full numerical linear stability results. The 
agreement in this case is very good (there are no degeneracies and associated higher-order contri- 
butions that may deteriorate the quality of the agreement as in the vortex breather case above); 
in fact, in some of the cases, the curves are almost indistinguishable throughout the considered 
parameter range. 

An important observation concerns, however, the role of the "hard" nature of the potential. 
In particular, as illustrated in Fig. 4, the quantity du/dJ is positive in this case, i.e., its sign 
is opposite from the soft case of the Morse potential (where doo/dJ = — 1). This results in the 
corresponding reversal of the stability conclusions in Figs. 5 and 6, in comparison with Figs. 1 and 
2 of the Morse case. That is, in-phase modes are now stable, while out-of-phase ones are unstable 
(as is true for the defocusing nonlinearity DNLS case also), while the reverse was true in the Morse 
potential (as well as for the focusing DNLS case). Lastly, we recall that since this is an even 
potential and thus f(ir/2) = 0, the leading order calculation would yield a vanishing contribution 
to the eigenvalues for the vortex case and a higher-order calculation is necessary to determine the 
stability of the latter. 

As an aside towards obtaining a fully analytical prediction for this case (as some of the quantities 
need to be obtained numerically above), we note the following. Although we cannot acquire an 
exact form for J (to), as in the case of the Morse potential, an approximate form for J can be found 
by using the rotating wave approximation (RWA), i.e. by supposing that x(t) « 2z% cos(wt). The 
introduction in the dynamical equations for the single oscillator leads to: 




(57) 




(58) 
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Figure 5: (Top panels) Profiles of an in-phase (left) and an out-of-phase (right) 2-site breather 
with the hard (j) 4 potential; u) = 3 and e = 0.05. The bottom panels show the dependence of the 
characteristic exponents a, of the corresponding configurations, on the coupling parameter e. The 
dashed lines correspond to the predictions of the stability theorems and dash-dotted lines to the 
RWA predictions, while the solid ones represent the full numerical result. 



UJ 2 



^i = Y^— (59) 

Thus, J = 2u)z 2 = 2uj(uj 2 — l)/3 and<9u;/<9J = 1/[2(cj 2 — 1/3)], and the corresponding expressions 
for the eigenvalues read: 

(7(0) « ^-e^f-^icos&xM (60) 



uj 2 — 1 

o-(tt) ~ ]l e^—j{cos(p)x q (0) (61) 

A comparison between the numerically acquired values of J{oo) and doo/dJ with the ones calcu- 
lated from the RWA is shown in figure 4. The agreement is remarkable and attests to the quality of 
the "single frequency" rotating wave approximation. In Figs. 5 and 6 the characteristic exponents 
calculated numerically (solid lines) as well as using Eqs. (57)-(58) (dashed lines) and via Eqs. 
(60)-(61) are compared, illustrating the excellent agreement between all three. 
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Figure 6: Same as Fig. 5, but for the three-site in-phase (left) and out-of-phase (right) configuration. 
Again the dash-dotted lines in the bottom panels represent the (fully-analytical) RWA predictions, 
which agree well with the semi-analytical dependences (dashed lines) and, in turn, with the full 
numerical results (solid lines). 



5 Conclusions and Perspectives 

The results presented in this work underscore the formulation of a toolbox that enables the sys- 
tematic characterization of both the qualitative and even the quantitative aspects of stability of 
multibreather and vortex breather waveforms in these large number of degree of freedom, Hamil- 
tonian lattice systems of the Klein-Gordon variety. A systematic calculation of the corresponding 
Floquet multipliers is presented and highlights the crucial components that imply stability, namely 
the proper combination of the sign of the coupling constant, the nature (hard or soft) of the poten- 
tial and the relative phases between the adjacent excited sites. E.g., for positive couplings, and soft 
potentials, out-of-phase structures may be stable near the vanishing coupling limit, while in-phase 
ones are unstable; the nature of the conclusions is reversed for either (small) negative couplings or 
for hard potentials. The explicit analytical predictions have been tested against numerical results 
both for symmetric (such as the hard <f) ) and asymmetric (such as the Morse) potentials, both for 
hard and soft ones, and both for simpler, non-degenerate one-dimensional multibreather settings 
and for more complex and degenerate two-dimensional vortex breathers. In all cases, the two the- 
ories whose results were shown to be equivalent herein, namely the Aubry band theory and the 
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MacKay Effective Hamiitonian method yield excellent qualitative and good quantitative agreement 
with the full numerical linear stability results. The latter may not be true only in degenerate cases 
where higher order contributions may be critical in breaking the relevant degeneracy (as we saw in 
the case of the discrete vortices for the Morse model). 

Naturally, a number of interesting directions for future consideration hereby arise. Perhaps 
the canonical one among them would involve a systematic derivation of higher order corrections 
for prototypical cases where the leading order approach yields vanishing results. For instance, 
the characterization of the stability of discrete vortices in the "super-symmetric" case of phase 
difference <f> = ir/2 for even potentials would be a natural example. Another possibility that is also 
emerging and would be relevant to consider from a mathematical point of view would be to examine 
models with inter-site nonlinearities, such as ones of the Fermi-Pasta-Ulam type. In these cases, 
where the potential is a function V(x n — x n -i)> it is relevant to point out that upon consideration 
of the so-called strain variables r n = x n — x n -i, the problem is reverted to an on-site potential case, 
for which it would be worthwhile to explore methods similar to the ones analyzed herein. These 
directions are presently under consideration and will be reported in future publications. 
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